Accuracy of BCS-based approximations for pairing in small Fermi systems 
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We analyze the accuracy of BCS-based approximations for calculating correlation energies and 
odd-even energy differences in 2-component fermionic systems with a small number of pairs. The 
analysis is focused on comparing BCS and projected BCS treatments with the exact solution of the 
pairing Hamiltonian, considering parameter ranges appropriate for nuclear pairing energies. We find 
that the projected BCS is quite accurate over the entire range of coupling strengths in spaces of 
up to about ~ 20 doubly degenerate orbitals. It is also quite accurate for two cases we considered 
with a more realistic Hamiltonian, representing the nuclei around 117 Sn and 207 Pb. However, the 
projected BCS significantly underestimates the energies for much larger spaces when the pairing is 
weak. 



I. INTRODUCTION 

A theory of nuclear pairing based on the BCS ap- 
proximation was considered for the first time 50 years 
ago [lj. Since then, BCS-based approximations or more 
generally the Hartree-Fock-Bogoliubov equation have be- 
come ubiquitous for calculating nuclear energies in the 
framework of density functional theory or self-consistent 
mean-field theory. It is thus important to understand the 
limitations of these approximations and use more accu- 
rate theory when needed. In particular, the BCS Ansatz 
of a condensate with indefinite particle number becomes 
problematic in finite systems with weak pairing, as is 
the case for nuclei at shell closures. Several ways have 
been proposed to improve the theory 0, beginning with 
number-projected BCS (PBCS) first proposed by Bay- 
man Q and Blatt Q. It is our aim here to evaluate 
the PBCS by testing it in situations for which an exact 
solution is available. 

To investigate the accuracy of PBCS approximation we 
shall consider a finite number of spin-1/2 fermions, e.g., 
neutrons or protons, distributed in a sequence of single- 
particle levels and interacting through a pairing force. 
We will mainly consider the reduced BCS Hamiltonian 
given by 

n n 
H = ^£i{a\ai + aUj) - s^ajotaja,-. (1) 

i i,3 

Here g is the strength of the pairing force acting in a 
space of f2 two-fold degenerate orbitals with the single- 
particle energies £j. 

The Hamiltonian (1) is exactly solvable 0, [|| and it 
was used in 60ies to make a critical analyses of the BCS 
approximation in finite Fermi systems. Thus in Ref. 6] 
Richardson studied the exact and the BCS solutions of 
the Hamiltonian (1) with = i and for systems with 
£1 = 8 — 32 at half filling, i.e., with the number of parti- 
cles N = fi. Such systems plausibly model the pairing in 
deformed nuclei with the active nucleons (for pairing cal- 
culations) filling the major shells 8-20, 50-82 and 82-126. 
The main conclusion of Ref. [(J was that BCS model 



strongly underestimate the pairing correlations even for 
relatively large values of the pairing strength. The ques- 
tion we address in this study is how much one could im- 
prove the BCS results relative to the exact model if we 
perform PBCS calculations. Some of the issues analyzed 
in this paper where also discussed recently in relation to 
metallic grains studies [j| [t| [l(| • There has also been a 
recent study in the nuclear physics context [TTI | . These 
authors found a significant difference between the exact 
solution and the PBCS. In Sec. II we shall argue that the 
PBCS approximation is nevertheless quite accurate if it 
is applied in a limited window around the Fermi level of 
the order of one major shell in atomic nuclei. 

Unfortunately, Richardson's model requires that the 
interaction matrix elements be equal in the pairing 
Hamiltonian It is essential to be able to treat the 

most general form of the Hamiltonian (1), with the ma- 
trix elements computed as integrals over an effective two- 
particle potential, if the theory is to be a global one de- 
scribing the entire nuclear mass table. In that case the 
Hamiltonian has a more general form 

si n 

h = Y^ £i(4 a i + °\ a i) - v ij a l a l a ] a 3> ( 2 ) 

i i>j 

where Vij = v^jj are calculated with some two-body in- 
teraction such as in Eq.(I6) below. There is no algebraic 
solution for this more general case, but we can obtain use- 
ful results up to and beyond f2 = 16 using ordinary nu- 
merical matrix methods. Some examples will be treated 
in Sec. IIIIl The realistic calculations show that the 
PBCS approximation gives accurate results, confirming 
the conclusions based on Hamiltonian (I) and Richardson 
model. The fact that PBCS can provide a good descrip- 
tion for the ground state of realistic Hamiltonians can be 
also seen from the large overlaps betweent the PBCS and 
the shell model wave functions [HI]. 

A good agreement between PBCS and the exact 
Richardson's solution we also find for the occupation 
probabilities of single-particle levels, analysed in Sec. IV. 
It is worth to emphasize that this agreement is obtained 
with two different wave functions for the ground state 
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of the systems, i.e., a condensate formed by identical 
Cooper pairs in PBCS model, a non-condensate structure 
based on distinct pairs in the case of exact Richardson's 
solution. These differences manifest clearly in the col- 
lectivity of the Cooper pairs, discussed in Sec. IV. The 
reason why such differences do not affect significantly 
the correlation energies and the occupation probabilities 
when the calculations are done in limited window around 
the Fermi level is not yet clear [l4j ■ 

Finally we would like to mention that another method- 
ology that is widely used to go beyond the BCS theory is 
the Lipkin-Nogami approximation. We do not consider it 
here for following reasons. First, it has been thoroughly 
studied in the past and its strengths and deficiencies are 
well known. It has a serious shortcoming in that the ap- 
proximation is not reliable near closed shell nuclei[16j|. 
Since we seek approximation methods that cover all the 
extremes that arise in the nuclear mass table, we find this 
method unsuitable. 



II. SOLUTIONS OF THE BCS HAMILTONIAN 



As shown in Ref. the pairing Hamiltonian (1) can 
be solved exactly. The solution resembles eq. ([5]) except 
that the operator is replaced by N pa i r different pair 
operators B u , 



/V 



i*>=n^i°)- 



The pair operators have the form 

1 
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They depend on energy parameters E v obtained by solv- 
ing the set of nonlinear equations 
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= 0. 



(8) 



The sum of pair parameters E v gives the total energy of 
the system, i.e., 



(9) 



The Hamiltonian (1) certainly describes the basic fea- 
tures of nuclear pairing correlations, and it is commonly 
solved using the BCS or PBCS approximation. Both 
BCS and PBCS methods are variational in that the ap- 
proximation is made on the wave function, and the en- 
ergy is calculated as an expectation value. The usual 
form of the BCS wave function is the well-known expres- 
sion np(iii+«iotat)|), but for putting it in the context of 
the other treatments one can write it as an exponentiated 
product of a pair operator, 
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The BCS ground state can then be expressed as a coher- 
ent superposition of pairs, i.e., 

|5CS)oc e r >> = £^|0>. (4) 

n 

The mixing amplitudes of the pair operator, written usu- 
ally as Xi — Vi/u.i, are given by the well-known BCS 
equations. The PBCS approximation is obtained by re- 
stricting the expansion in Eq.(3) to the term having the 
required number of particles. Thus, in PBCS the ground 
state wave function can be expressed 



\pbcs) oc (r^—io), 



(5) 



where N pa i r is the number of pairs. The PBCS equations, 
which determine the mixing amplitudes x% of the pair 
operator (2), are derived by minimizing the average of 
the Hamiltonian in the state (5). They can be solved by 
using the residual integrals technique [151 ] . Alternatively, 
if the number of pairs is not too large, the amplitudes 
Xi can be found by using the technique of recurrence 
relations. 



In the limit g = the pair energies E v of the ground 
state solution coincide with the lowest single-particle en- 
ergies, i.e., E v = 2e„, (y = 1, 2, ...N pair ). When the inter- 
action is turned on, the pair energies evolve toward lower 
values and could become complex two at a time. This 
fact was used by Richardson to obtain a set of equations 
in which the singularities are removed [f|. 

For small values of g in finite systems, the BCS con- 
densate collapses and the BCS approximation gives zero 
correlation energy. On the other hand, due to the fi- 
nite distance between the levels, for small values of the 
interaction strength the pairing interaction energy can 
be calculated perturbatively. This is opposite to what 
happens in finite systems where the pairing correlations 
depends exponentially on pairing strength in the weak 
coupling limit. The perturbative solution can be easily 
derived by enumerating the two-particle two-hole config- 
urations starting from the lowest energy configuration, or 
by taking the small-g limit of the Richardson equations 



The second-order perturbation result for the interac- 
tion energy, i.e., the energy gained by the system when 
the interaction is turned on, is given by 
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This expression is valid for even N; for odd N = N pa i r + 1 
the j sum begins at j = N pa i r + 2. It can be shown Q 
that the approximation (|10[) for the interaction energy is 
valid if g < gp = g*(l — g*), where g* is the convergence 
radius of the perturbative expansion given by 
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FIG. 1: Correlation energies for the Hamiltonian (1), calcu- 
lated in various approximations. Upper panel: Q — 8, N pa i r — 
4; lower panel: SI = 80, N pair = 40. 



FIG. 3: Errors for the correlation energies calculated in the 
PBCS approximation. 




A. Correlation Energies and Pairing Gaps 

In this section we shall discuss the correlations energies 
and the pairing gaps provided by the Hamiltonian (1) for 
a single-particle spectrum formed by Q two-fold degen- 
erate orbits with uniformly spaced energies, i.e., £< = i. 
The analyses is done for systems formed by N particles 
distributed in Q = N levels with O = 8 - 80. For the 
strength of the interaction we use the range g = 0.1 — 0.8, 
which covers all the interesting coupling regimes met in 
nonspherical atomic nuclei. In all figures discussed in 
this section the energies, the pairing gaps and the inter- 
action strength are given in units of single-particle levels 
spacing. 

We first discuss the correlation energies calculated ex- 



actly and in various approximations. The correlation en- 
ergies are defined by the difference in energies between a 
single Slater determinant and the pair-correlated state, 
i.e., 

E corr (g) = E HF -E(g). (12) 

The results are shown in Fig.l. From this figure is obvi- 
ous that the BCS seriously underpredicts the correlation 
energy while the PBCS does much better. To see the dif- 
ferences more quantitatively, in the next figures we show 
for each approximation how the error depends on g and 
f2. The BCS errors shown in Fig. 2 are large, making this 
approximation completely unreliable. The results for the 
PBCS approximation are shown in Fig. 3. The errors 
are much smaller, but become unacceptably large for the 
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biggest space, fl = 80. In spacing corresponding to a 
single major shell the error is well under 10%. It is im- 
portant to emphasize that the example discussed above 
does not contradict the fact that pairing models become 
more accurate if one moves toward the thermodynamical 
limit [7(. To reach this limit the particle number is in- 
creased but, in the same time, the calculations are done 
keeping a fixed energy window around the Fermi energy. 
In this way the increase of particle number has as ef- 
fect an increase of level density around the Fermi level 
which, in turn, is increasing the pairing correlations and 
by that the accuracy of pairing models. This is different 
from what happens in the system with 40 pairs discussed 
above and in many BCS and HFB calculations performed 
for atomic nuclei, where the energy window around the 
Fermi level is increased such that to include the deep 
bound nucleons. As discussed above, by this procedure 
the accuracy of pairing calculations becomes worse not 
better. 

The question which arises is why BCS approximation 
strongly underestimates the pairing correlations in finite 
systems. This question was recently addressed in rela- 
tion to metallic grains calculations. Thus, in Ref.[9( it is 
argued that BCS works properly only for so-called " con- 
densed" levels, i.e., the levels with the energies located in 
the interval I = |A — where A is the gap parameter 
in the BCS equations and fi is the Fermi energy. This 
conclusion is supported by the observation that the cor- 
relation energy calculation in BCS is close to the exact 
result obtained if from the exact solution is kept only the 
contribution of condensed levels (in the exact solution the 
condensed levels have usually complex pair energies) . On 
the other hand, it was found that the contribution of lev- 
els located outside the interval / are underestimated in 
BCS. Since in the exact solution the Cooper pairs cor- 
responding to these levels have the pair energies close 
to the single-particle levels, one expects that the con- 
tribution of these levels to pairing correlations could be 
treated perturbatively. Based on these observations it 
was found [1, [l(| that in metallic grains the correla- 
tion energy could be approximated by a formula com- 
bining the BCS and the perturbative expressions in the 
sum Ecorfig) + A + a(g)E^ orr , where the last term is the 
perturbative result and a(g) is a function of the order 
of unity determined by the fitting protocol. The BCS 
contribution is approximated by the first two terms rep- 
resenting the condensation energy and the pairing gap for 
infinite system. The later accounts for the size correction 
to the bulk result (first term) . In order to see if such an 
approximation could work for the small systems analyzed 
here, we have simply replaced the first two terms in the 
above expression by the BCS result (for finite systems) 
and we take a(g) = 1, 
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FIG. 5: Pairing gaps (Eq. (|14[) ) calculated in various ap- 
proximations. Upper panel: = 8, N pa ir = 4; lower panel: 
Q = 80, N pair = 40. 
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FIG. 6: Error in pairing gaps calculated with the PBCS 
approximation. 



works surprisingly well for a wide range of the pairing 
strengths, including the physical region up to g ~ 0.8. 
Of course, it is not applicable to situations where there 
is a degeneracy in the single energies at the Fermi level, 
since perturbation theory diverges in that situation. 

We next compare the pairing gaps calculated in the 
BCS and PBCS approximations with the exact values. 
The gap is defined as the second difference of energies 
between an odd-iV system and its neighbors with even- 
N, 



p _ pBCS , pP 



(13) 



The corresponding results are shown in Fig. 1 by in- 
verted triangles and the corresponding errors are given 
in Fig. 4. It can be seen that this simple approximation 



A (3) (N) = - (2E(N) - E(N - 1) - E(N + 1)) . (14) 

We shall call A^ 3 ^(iV) the gap at number N. The results 
at the smallest and largest spaces are shown in Fig. [5] 
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One can see that PBCS gives accurate gaps in all cou- 
pling regimes, contrary to the larger systems analyzed in 
metallic grain studies [1, H, The errors associated 
with the PBCS approximation are shown in more detail 
in Fig. [6] It is probably acceptable to tolerate an error 
up to 10% for the gap, but not higher. This confirms for 
the gaps that the PBCS is only reliable up to moderate 
size spaces. 

It is interesting to note that in spite of the large er- 
rors in the correlation energies, in the physical region 
of the strength parameter the BCS pairing gaps come 
much closer to the exact results. However, contrary to 
the PBCS results, the good agreement of the BCS gaps 
to the exact values in the region of well-developed pairing 
correlations is just a manifestation of the errors cancella- 
tion when the subtraction is performed in equation (14). 
Consequently, fixing the pairing force by the odd-even 
mass difference, as usually done in nuclear structure cal- 
culations, does not guarantee a good description of cor- 
relation energies within the BCS approximation. 



III. PBCS FOR REALISTIC PAIRING 
INTERACTIONS 

Ultimately, the theory of nuclear binding energies 
should be based on realistic interactions dropping the 
constant-g approximation. The Richardson solution has 
been somewhat generalized to encompass separable pair- 
ing interactions 26], but to be truly realistic the Hamil- 
tonian must allow completely general interactions in the 
many-body space of pairwise occupated orbitals. These 
Hamiltonians can be solved by straightforward configu- 
ration interaction (CI) methodology, in which a Hamil- 
tonian matrix is constructed in the Fock space of the 
orbitals and diagonalized by standard linear algebra op- 
erations. The size of the space D needed to represent the 
most general paired wave function is given by the number 
of combinations of N pa i r orbitals out of total of f2, 



D = 



(15) 



For example, for fl = 16 orbitals and N = 16 particles 
the dimension of the space is 12,870. The lowest eigen- 
vector for a space of this size is easily calculated on serial 
computers using the Lanczos algorithm. 

We have carried this out for two examples in which the 
pairing is very different. The Hamiltonian makes use of 
the orbital energies and wave functions from the global 
mean- field calculations of Ref . [17J , which are based on 
the Skyrme SLy4 energy functional. Pairing is strong 
in the Sn isotope chain, and we will take g^Snso and 
its neighbors as an example of where the pairing is well 
developed. The second example is 207 Pb. In the global 
systematics of neutron pairing gaps, the one at 207 Pb is 
the smallest (A' 3 ' = 0.32 MeV), so this provides a good 
test of the approximation methods in the weak pairing 
limit. 



We use the published code ev8 [18| to recalculate the 
needed orbital properties, starting from the wave func- 
tions provided in the original global survey [3. The 
orbitals are represented internally in the code with a 3- 
dimensional mesh, so they need not have good angular 
momentum quantum numbers. However, for the even-iV 
Sn and Pb isotopes, the mean field solution is spherically 
symmetric and the orbitals can be given angular momen- 
tum assignments. In Ref. [17| the pairing interaction 
was taken as a density-dependent contact interaction in 
a space truncated to a band of width 10 MeV about the 
Fermi energy. Here we use an ordinary delta function to 
generate the pairing matrix elements, 



Vij = v / (fV|(/>j(r)| 2 |<^-(r)| 



(16) 



where the cj>(r) are orbital wave functions. The strength 
Vo has been fitted to the global systematics of pairing 
gaps and the self-consistent orbitals were generated 
with that value. The single-particle energies and the ma- 
trix of values were then used as input data for separate 
codes to solve the Hamiltonian Eq. © . We calculate the 
correlation energies and the pairing gaps using a range 
of values for vq obtained by scaling the matrix elements 
obtained from the ev8 code. 

For the Sn isotopes with neutron number N around 
68, there are fl — 16 orbitals in the 10 MeV window, 
originating from the d 5 / 2 , 97/2, S1/2, d 3 / 2 and h n / 2 shells 
of the spherical shell model. This space is small enough 
to permit the exact calculations to be performed with- 
out special numerical difficulties. Fig. [JJ shows the neu- 
tron correlation energy in 116 Sn as a function of inter- 
action strength vq- The range for vq includes the value 
vo ~ 450 MeV fm 3 that has been fitted to the global gap 
systematics using the BCS approximation [l9(. The com- 
parison in Fig. [7J confirms the results obtained with the 
Hamiltonian with the constant-g pairing. Namely, the 
BCS systematically underpredicts the correlation energy 
while the PBCS tracks the exact energy very well. PBCS 
describes also very well the pairing gaps, as shown in Fig. 
8. There is also a good agreement for the BCS gaps but, 
as already noticed in the previous section, this agreement 
is in fact a manifestation of errors cancelation. 

Next we show the results for Pb. Here the energy trun- 
cation to a 10 MeV window gives a space that is still too 
large to easily perform the exact calculation, so we trun- 
cated it further (~ 8 MeV window) to leave = 16 or- 
bitals. Fig. [5] shows the correlation energy in this space 
and the BCS and PBCS approximations. One sees that 
the PBCS keeps an accuracy much better than 100 keV, 
while the BCS is off by more than a half of an MeV. 
The correlation energies in the other nuclei needed for 
the 207 Pb gap are very small and in BCS there is no con- 
densate in 207 Pb and 208 Pb. As a result, the BCS error 
for the 206 Pb correlation energy is not well canceled in 
formula for the gap energy. Fig. [TTJ] shows the calculated 
gap in the three treatments. One sees again that the 
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FIG. 7: Pairing correlation energy in 116 Sn. The pairing in- 
teraction is a delta function with a strength scaled from the 
nominal value Vo = 465 MeV-fm 3 by a factor given on the ab- 
scissa. The orbital space is the 16 orbitals around the Fermi 
energy as described in the text. Solid line shows the result of 
exact diagonalization. 



FIG. 9: Pairing correlation energy in 206 Pb. For these calcu- 
lations, the space was truncated to Q = 16 by including the 
11 highest orbitals below the N = 126 magic number and the 
5 gg/2 orbitals above N = 126. This corresponds to an energy 
window of about 8 MeV . 





v/v„ 



FIG. 8: Neutron pairing gap at 117 Sn. The correlations ener- 
gies for the three nuclei needed for Eq. (|14|l were calculated 
with the same functional and in the same space as in Fig. [7] 



FIG. 10: Neutron pairing gap at Pb. For these calculations 
we used the same space as in Fig. 9. 



PBCS is remarkably accurate. The BCS error of 100-200 
keV is quite significant on the scale of the empirical gap 
energies, which fluctuate around an average of 1 MeV 
with an r.m.s. deviation of 300 keV. 



the pairing force. The results discussed in this section 
correspond to the Hamiltonian (1) with Si = i. 



A. Occupation probabilities 



IV. OCCUPATION PROBABILITIES AND 
TWO-BODY CORRELATIONS 

To probe the accuracy of BCS-based models relative to 
the exact solution we have also analyzed the occupation 
probabilities and the two-body correlations induced by 



The BCS occupation probabilities are calculated by 
solving the standard BCS equations (here we take into 
account the renormalisation of the single-particle energies 
by the diagonal term of the interaction) while the PBCS 
values are obtained by using the residual integral method 
described in Ref. [la]. 

In the exact solution the occupation probabilities vf 
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FIG. 11: Occupation factors «f = vfil — vf) for N pair — 8, 
£1 — 16 and p = 0.32. The results correspond to the Hamil- 
tonian (1) with Ej = i. 



FIG. 12: The same as in Fig. 11 but for g = 0.42 



are given by 



N 

V 



^ (2i - E v f 



(17) 



where a v are obtained by solving the set of equations 
1 



1 (18) 



and C v are given by 



C u = 



2N 

E 



i 



(2» - 2?„ 



(19) 



To solve the equations above is convenient to rewrite 
them in terms of the real and the imaginary parts of 
pair energies E v . The corresponding expressions can be 
found in Ref. [6J. 

For the discussions below we shall use the product be- 
tween the occupation and non-occupation probabilities, 
i.e., nf = vf(l — vf), which provides relevant informations 
on the diffusivity of the Fermi sea and the entanglement 
properties of pairing tensor (see Eq.(21) below). 

In Figs. 11-13 are shown the values of /tf for a system 
with N=8 pairs and for three values of the strength pa- 
rameter corresponding to the weak (g = 0.32), interme- 
diate (g = 0.42) and strong coupling (g = 0.87) regimes. 
For these strength values the gaps (exact results) are ap- 
proximatively equal to 0.6, 1.0, 5.0, respectively. As seen 
from Figs. 11-13, the values of k| for the states outside 
the interval I = |A — fi\ are rather well described by 
PBCS and underestimated by BCS. On the other hand, 
one can notice that for the states which are the clos- 
est to the chemical potential BCS gives in the weak and 
intermediate coupling regimes larger values than PBCS 




FIG. 13: The same as in Fig. 11 but for g = 0.87 



and the exact solution. In the strong coupling regime 
shown in Fig. 13, when all levels are inside the inter- 
val /, the BCS and PBCS results become close to the 
exact values for all levels included in the calculations. 
Hence, when the pairing correlations are well-developed, 
BCS describes rather well the occupation probabilities 
even though the errors for the condensation energies are 
large. 



B. Two-body correlations 

The distribution of occupation probabilities around the 
Fermi level determines the intensity of two-body correla- 
tions, which could be eventually probed in pair transfer 
processes. Two-body correlations are commonly intro- 
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duced by the two-body density denned by 

p 2 {xi,x 2 ) = ^2 J |*(zi,Z2, ■ ■■,x N )\ 2 dr 3 ...dr N 
= (Q\a + (x 1 )a + (x 2 )a(x 2 )a(x 1 )\Q), 

where VP is the many-body wave function, x denotes the 
radial and spin coordinate, i.e., x = ra, and af(x) is the 
particle creation operator. To separate the genuine two- 
body correlations induced by pairing correlations from 
the correlations of Hartree-Fock type, the two-body den- 
sity is usually written in the following form [2Cl | 

p 2 (xi,x 2 ) = p(xi)p(:c2)-|p(:ri,X2)| 2 +|K(xi,:z;2)| 2 , (20) 

where p{x) is the (local) particle density while p(xi, x 2 ) is 
the non-local (exchange) part of particle density. The last 
term defines the genuine two-body correlations, i.e., the 
correlations not included in the independent mean-field 
picture of fermion motion. In the BCS approximation 
the last term corresponds to the pairing tensor in the 
coordinate representation, i.e., 

k{xi,x 2 ) = (0|o(a;2)a(a!i)|0) = Kiipi(xi)(pi(x 2 ), 

(21) 

where /t, = (0|aj<ii|0) = UiVi is the pairing tensor in a 
single-particle basis defined by the operators a] and ifii 
are the associated wave functions. 

According to its definition, the pairing tensor in the 
coordinate representation provides information about the 
spatial correlations between two nuclcons irrespective if 
they belong or not to the same Cooper pair. The spatial 
structure of these correlations in atomic nuclei have been 
recently discussed in Refs. [HI, H3|- If we need to inves- 
tigate only the spatial correlations between the nucleons 
belonging to the same pair, instead of pairing tensor one 
should analyse the pair wave function. The latter has a 
different structure in BCS-based models compared to the 
exact solution. Thus, in BCS and PBCS models all pairs 
are described by the same wave function 

4>(xi,x 2 ) = Cy j Xiipi{xi)ipj{x2), (22) 

i 

where the mixing amplitudes Xi and the normalization 
factor C depend on the approximation used to describe 
the condensate. Thus, in BCS Xi — Vi/ui and C = 
^^/w 2 , where Vi,Ui are the occupation amplitudes 
provided by the BCS equations. Formally, the same ex- 
pressions can be used for the PBCS model but in this case 
the amplitudes Vi,Ui are just variational parameters, not 
occupation amplitudes. 

For the exact solution each pair is described by a spe- 
cific wave function, i.e., 

<ft„(si,s 2 ) = 2s- - E ip i( x i) < Pi( x 2)- ( 23 ) 
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FIG. 14: The Schmidt numbers (Eq. 24) for a system with 
N pa i r = 8 and Q — 16. 



Therefore the correlations among the nucleons belonging 
to the same pair can be rather different according to the 
pair function they belong to. To characterize the amount 
of correlations in the pairing tensor (21) and in the pair 
wave functions (22-23) we employ the quantity [23| 

*T = l/£>f, (24) 

i 

where Wi are the mixing amplitudes of the normalized 
two-body wave functions (21-23). Th e q uantity K is 
sometimes called the Schmidt number [2J] and it gives 
a global indication of the degree of entanglement in two- 
body systems. Thus, the minimum value of K is 1 and 
is obtained when in the expansion (24) only one term is 
non-zero; this case corresponds to no entanglement since 
the two-body wave function is split in the individual wave 
functions of the two- particle system. The maximum pos- 
sible value of K is obtained when the weights Wi have the 
same value for all terms in the expansion. 

In Fig 14 is shown how K number evolves as a function 
of the strength parameter. One can see that the entan- 
glement properties of Cooper wave functions are rather 
similar in BCS and PBCS approximations and very dif- 
ferent from the entanglement of the correlation function 
(21). As already stressed above, the latter describes the 
correlations between two generic fermions, which do not 
necessarily belong to the same Cooper pair. 

A particular behavior have the Cooper pairs (23) de- 
scribed by the exact solution. Their entanglement prop- 
erties depends on how far are their pair energies E v (more 
precisely, their real part) from Fermi level. Thus, as seen 
in Fig 14, the 8th pair (i.e., the one corresponding to E$), 
which is the closest to Fermi level, is the most entangled. 
At the other extreme is the 1st pair, corresponding to 
Ei, which remains almost uncorrelated in all coupling 
regimes. An intermediate behavior has the 5th pair which 
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remains uncorrelated up to g w 0.53 and then its entan- 
glement is increasing rather fast, in a similar fashion as 
for the 8th pair. The strength value for which K num- 
ber of the 5th pair starts to increase corresponds to the 
value at which the pair energy E$ becomes complex. As 
mentioned previously, the pair energies become complex 
when the corresponding level (in the limit g = 0) becomes 
condensed, i.e., enters in the interval I = |A — fi\. Hence, 
in the exact solution at a given value of the strength 
only some pairs become entangled, namely the ones cor- 
responding to the condensed levels. This is very different 
from what happens in BCS and PBCS models in which 
all Cooper pairs are identical and therefore all of them 
have the same entanglement properties. 

V. SUMMARY AND CONCLUSIONS 

In this paper we have studied how reliable are the ap- 
proximations used to solve the pairing Hamiltonian when 
the parameters are chosen in a range appropriate to cal- 
culate nuclear ground states. As is well known, the BCS 
approximation does not do well for correlation energies. 
On the hand, we find that PBCS, with its variation after 
number projection, is highly accurate over most the inter- 
esting parameter range for moderate size orbital spaces. 
This applies to orbital spaces such as a single major shell 
or energies truncations of the order of 5 MeV around the 
Fermi level. Much larger spaces, for example including 
all the occupied orbitals, give a degradation in the per- 
formance of the PBCS that we do not yet understand. 
However, by renormalizing the pairing interaction, the 
calculations in large spaces could be reduced to smaller 
ones in which the PBCS gives better results. One way 
the renormalization could be done is to demand the same 
computed gaps when the number of states available for 
the active pairs is changed. For example, as seen in Fig. 
3, if in the PBCS calculations we decrease the number 
of active pairs from 40 to 8 and, in order to keep the 
same gap, we increase the strength from 0.27 to 0.5, the 
error in the correlation energy drops from 20% to below 
5%. Hence, to reduce the errors of PBCS calculations 
is preferable to restrict the calculations to a small num- 
ber of active particles with energies located close to the 



Fermi level. In any case, single-major shell truncations 
are widely used and there is no reason to not use the 
PBCS with those conditions. 

Both the BCS and PBCS seem to do well on calculating 
pairing gaps with the the reduced BCS Hamiltonian Eq. 
(1). However, in the case of BCS the improvement arises 
from a cancellation of errors. The BCS underpredicts 
the correlation energies, but more seriously is missing 
the true values for odd systems. In the realistic case of 

Pb the errors did not cancel well, and so we regard 
the BCS as unreliable at the level of 100 keV or so. The 
PBCS maintained its accuracy for the two realistic cases 
we considered, confirming our overall assessment of its 
reliability. 

We have also found that PBCS describes accurately 
the occupation probabilities of the single-particle lev- 
els. However, contrary to the agreement found for the 
correlations energies and the occupation probabilities, 
the entanglement properties of pair wave functions are 
very different in PBCS and in the exact solution. Thus, 
while in PBCS the entanglement of Cooper pairs depends 
smoothly on the interaction strength, in the exact wave 
function, formed by non-identical pairs, the entanglement 
properties of Cooper pairs depend strongly on the posi- 
tion of their pair energies relative to chemical potential. 

Finally, we would like to mention that a ground state 
based on non-identical pairs, specific to the exact solu- 
tion, could be more appropriate for the description of 
loosely bound systems such as nuclei close to the drip 
lines. For such nuclei one expects that the properties 
of Cooper pairs formed by the valence nuclcons moving 
in loosely bound and continuum single-particle state to 
be rather different from the pairs formed by the deeper 
bound nucleons. 
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